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There were some requests about DERIVE ’s compatibility with WINDOWS 8. 1 can appease 
all WIN 8 users - and all of you who intend to change to WIN 8. There are no problems, 
Josef. 


DERIVE and WINDOWS 8 

Our member Gunter Schodl provided some information concerning DERIVE and 
WINDOWS 8: 

Hello Josef! 

Derive can be installed under Win 8 without any problems. Like under Win 7 the Help file 
cannot be accessed without a patch. 

What you also can do is using Hyper-V under Win8 (it must be installed as a Windows com- 
ponent using appwiz.cpl, activate hardware virtualisation in the Bios of the PC), then you can 
install a virtual WinXP or Win 7 and run DERIVE. 

The link for the patch (WinHlp32.exe) is 

http://www.microsoft.com/de-de/download/details.aspx?id=35449 

Then DERIVE will run as usual. 

Greetings 

Gunter 

Another valuable note from Gunter: 

There is a nasty Macro-error message appearing when calling the DERIVE Online-Help. You 
can find advice how to avoid this message at: 

http://support.microsoft.com/kb/9176Q7/de (German) 

http://support.microsoft.com/kb/917607/en (English) 

Latest news from TIME 2014 

It is a great pleasure to announce the list of keynote speakers for TIME 2014: 

Peter Baumgartner (Danube University Krems) 

Regina Bruder (Technical University Darmstadt) 

Bruno Buchberger (RISC Institute Linz) 

Pavel Pech (University of South Bohemia, Budejovice) 

✓ 

Gilles Picard (Ecole de technologie superieu re, Montreal) 

Marlene Torres -Skoumal (VIS and International Baccalaureate Schools) 
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Liebe bUG-Mitglieder, 

Ich musste (und wollte) bis heute auf ei- 
nen schon lange versprochenen Beitrag war- 
ten. baher die Versptitung mit bNL#88. 

Es war meine feste Absicht, die letzte 
Folge von Piotr Trebisz' Schneckenhausern zu 
bringen und war auch schon ziemlich fertig 
mit der Ubersetzung als Statistik 4 von Guido 
Herweyers eine ungeheure Eigendynamik be- 
kam. Einerseits konnte ich nicht widerstehen, 
die Simulationen mit bERIVE nachzuvollzie- 
hen und andererseits - wichtiger - gab es fur 
mich noch offene Fragen, die Guido in dan- 
kenswerter Weise sehr umfassend beantwor- 
tete. ber Anhang zu den Funktionen von Zu- 
fallsvariablen ist fur mich sehr wertvoll. 

Fred Tydeman hat im letzten bNL be- 
klagt, dass sich zu seinem in den bERIVE 
News vorgestellten Problem mit dieser spe- 
ziellen Folge kein bERIVIAN gemeldet hat. 
Nun, in diesem bNL gibt es gleich drei z.T. 
sehr ausfiihrliche Beitrage zu dieser Folge. 

Es ist erfreulich, dass gelungen ist, Jo- 
hann Wiesenbauer zu einer neuen Ausgabe 
seine Titbits zu ermuntern. In einem zweitei- 
ligen Aufsatz beschreibt er eine bERIVE 
Implementierung des nach dem RSA- 
Algorithmus bekanntesten Public-Key Ver- 
schlusselungssystems. 

Bitte beachten Sie die wertvollen Hinwei- 
se zu bERIVE & WIN 8, die uns Gunter 
Schodl zur Verfugung stellt. Er ist da immer 
sehr rasch voll informiert - und damit wir mit 
ihm, herzlichen bank dafiir lieber Gunter. 

Freuen Sie sich mit mir auf die ntichsten 
Ausgaben mit u.a. Beitrtigen uber Primzahlen 
und Taylorreihen (b. Oertel), das Briisseler 
Tor (E. van Lantschoot) und naturlich liber die 
Schneckenhtiuser, von denen Sie eine Kost- 
probe auf Seite 53 sehen konnen. 

Es bleibt mir noch, Ihnen alien ein gluckli- 
ches, erfolgreiches und gesundes Jahr 2013 
zu wunschen. 


bear bUG Members, 

I had (and I wanted) to wait for a long 
promised contribution until today. And this is 
the reason for the delay of bNL#88. 

It was my intention to include the last 
part of Piotr Trebisz' snail house series and 
its translation was almost ready. But then 
Statistics 4 (Guido Herweyers) got an im- 
mense self dynamic. I could not resist repro- 
ducing the simulations using bERIVE and at 
the other hand - more important - there were 
some open questions for me which were kindly 
answered very comprehensive. In my opinion 
the appendix to the functions of random vari- 
ables is very valuable. 

In the last bNL Fred Tydeman complained 
that no bERIVIAN responded to his special 
sequence problem presented in the bERIVE 
News Group. Now, in this bNL he - and all our 
members - can find three very detailed con- 
tributions. 

I am very happy that I was able to en- 
courage Johann Wiesenbauer to a new Tit- 
bits-contribution. He describes in two parts a 
bERIVE implementation of the famous EIGa- 
mal encryption algorithm which is together 
with the RSA-encryption the best known pub- 
lic-key encryption method. 

Please notify the valuable notes to bE- 
RIVE & WIN 8 which are provided by Gunter 
Schodl. He is always fully informed about 
hard- and software news. Many thanks for 
your support, dear Gunter. 

Be looking forward to the next issues 
with among others contributions on prime 
numbers (b. Oertel), the “Brussel's Gate" (E. 
van Lantschoot) and on the snail shells. You 
can have a tasting on page 53. 

For me remains wishing you a happy, suc- 
cessful and healthy year 2013. 

Viele GriiRe, kindest regards 


Download all Z)7 VL -Derive- and Tl-files from 

http : //www . austromath . at/dug/ 
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The DERIVE-NEWSLETTER is the Bulletin 
of the DERIVE & CAS-77 User Group. It is 
published at least four times a year with a 
content of 40 pages minimum. The goals of 
the DNL are to enable the exchange of ex- 
periences made with DERIVE, 77-CAS and 
other CAS as well to create a group to dis- 
cuss the possibilities of new methodical and 
didactical manners in teaching mathematics. 


Editor: Mag. Josef Bohm 
DXust 1, A-3042 Wurmla, Austria 
Phone: ++43-(0)660 3136365 
e-mail: nojo.boehm@pgv.at 


Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the DNL. It 
must be said, though, that non-English arti- 
cles will be warmly welcomed nonetheless. 
Your contributions will be edited but not 
assessed. By submitting articles the author 
gives his consent for reprinting it in the 
DNL. The more contributions you will send, 
the more lively and richer in contents the 
DERIVE & CAS -TI Newsletter will be. 

Next issue: March 2013 


Preview: Contributions waiting to be published 

Some simulations of Random Experiments, J. Bohm, AUT, Lorenz Kopp, GER 
Wonderful World of Pedal Curves, J. Bohm, AUT 
Tools for 3D-Problems, P. Luke-Rosendahl, GER 
Hill-Encryption, J. Bohm, AUT 

Simulating a Graphing Calculator in DERIVE , J. Bohm, AUT 

Do you know this? Cabri & CAS on PC and Handheld, W. Wegscheider, AUT 

An Interesting Problem with a Triangle, Steiner Point, P. Luke-Rosendahl, GER 

Graphics World, Currency Change, P. Charland, CAN 

Cubics, Quartics - Interesting features, T. Koller & J. Bohm, AUT 

Logos of Companies as an Inspiration for Math Teaching 

Exciting Surfaces in the FAZ / Pierre Charland 's Graphics Gallery 

BooleanPlots.mth, P. Schofield, UK 

Old traditional examples for a CAS - what’s new? J. Bohm, AUT 

Truth Tables on the TI, M. R. Phillips, USA 

Where oh Where is It? (GPS with CAS), C. & P. Leinbach, USA 

Embroidery Patterns, H. Ludwig, GER 

Mandelbrot and Newton with DERIVE , Roman Hasek, CZK 

Some Projects with Students, R. Schroder, GER 

Structures in the Set of Prime Numbers, D. Oertel, GER 

Dirac Algebra, Clifford Algebra, D. R. Lunsford, USA 

Laplace Transforms, ODEs and CAS, G. Picard & Ch. Trottier, CAN 

A New Approach to Taylor Series, D. Oertel, GER 

Cesar Multiplication, G. Schodl, AUT 

Henon & Co; Find your very own Strange Attractor, J. Bohm, AUT 

Rational Hooks, J. Lechner, AUT 

Mathematical Model for Snail Shells (4), P. Trebisz, GER 

Simulation of Dynamic Systems with various Tools, J. Bohm, AUT 

An APL-like SHAPE function in DERIVE 6, D. R. Lunsford, USA 

Brussels Gate in Dendermonde, E. van Lantschoot, GER 

Recursive Series of Numbers, Polynomials and Functions, D. Halprin, AUS 

and others 


Impressum: 

Medieninhaber: DERIVE User Group, A-3042 Wurmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag. Josef Bohm 
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Statistics with Tl-Nspire 3.1/3.2 (Part 4) 

Visualising and Simulating Dynamically with Tl-Nspire 

Guido Herweyers, KHBO Campus Oostende 
quido.herwevers@khbo.be 


Part 3: Discovering Probability Distributions 


(1) z-scores versus t-scores 


Set A' N(ju,a) and take a sample of size n from this population. Then the mean of this sample is given 

_ r 


by: X ~ N 


V 


4n , 


X — LI 

. Standardisation gives Z = — ~ 7V(0,1). 

<j 

4n 


In case of an unknown standard deviation o of the population o will be estimated by the standard de- 

X — u 

viation of the sample 5. Then we have a new random variable T = . 


x — LI x — LI 

We compare the “z-scores” z = — — — with the “t-scores” t = — - — using their density diagrams 

Jn 4n 

performing a simulation with the population mean ju= 100, the population standard deviation o = 10 
and a sample size ?/ = 4. (You can find hints how to reproduce the screen below in DNL#87 page 6.) 



-5-4-3-2-10123456 

zscores 


sample_size n.= 4 

sample ► { 88.3809, 104. 573,88. 699 1 , 1 15.026 } 
sample_mean:=meani sample I ► 99.17 
sample_mean- 100 


zscore: = 


-0.166005 


sample_mean- 100 

tscore:= : ► “0.127738 


stDevSampl. sample) 


dim^ zscores) ► 2012 


A sample 


B 


zscores 


c tscore* 


♦ =randnorm(1 00,1 0,4) =capture(zscore,1) =captur< 


1 88.3809 -1.78279 


2 104.573 0.970368 

3 88.6991 -0.325437 


4 115.026 -0.079782 

5 0.127264 


6 
7 

^[B5 =0.12726380326 


-0.764428 

1.16868 


The t-scores show a larger dispersion than the z-scores 



The t-scores are following a t-distribution with 3 degrees of freedom (3 = n- 1). Its graph is also bell 
shaped, but it is not the normal distribution! 


(Note: You can add the graphs of the distributions via Data & Statistics in the Documents Toolbox. 
Choose options 4: Analyze > 4: Plot Function.) 


(2) Functions of random variables 

Take any random number X in the interval [0, 
distributed following the density function 

/ 0 ) = 

If a random number between 0 and 1 is generated then the random variable X gets a certain value x 
(lower case!) according to the probability mechanism of the density function f 

Generating many random numbers in the interval [0, 1] and presenting the respective density histo- 
gram we should observe that this histogram will become more and more stable and finally “converge” 
to the density function (probabilities are relative frequencies for many tries: law of large numbers). 

200 simulations (200 random numbers generated) will show big variation in the density histogram:. 

In the following we show two simulations with 200 random numbers uniformly distributed in [0, 1], 
the next one is the presentation of 2500 random numbers and finally we will generate 5000 random 
numbers in the interval. 


1]. Then X (upper case!) is uniformly (or rectangular) 

Jl 0 < x < 1 
[0 else 
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sampSe size nt-300 
i -rand(n) 

* { ( J.07KT76 0,M |0ft5,0. ] 63ff7fi,C> J 77;36?.(li i S4»&59.0 < «»>|fi.0 4 566+57.0,4 1 1.0,017 1 1 1 ,ClB ] 5* 



For n = 2500 and then for n = 5000 variation becomes much smaller. The density histogram looks 
more and more like a rectangle (in case of not too small class widths). 




(Note: for n = 10000 we receive an error message: resource exhaustion.) 

In this way one can get a good impression of the probability distribution of a random variable by 
simulation. 

Variable X can be used to define new random variables, e.g. U = X 2 . The probabiliy distribution (the 
density function) of U is approximated by the density historgram of a large number of generated 
random values u = x 2 . 
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The density function of U is the first derivative of the distribution function F v of U: 
For 0 <x < 1: 

F v (x) = P(U < x) = P(X 2 <x) = P{-4x <X<yfx)= I bdt = yfx 

0 


Hence, the density function of U is given by: f v {x) = F^(x) = 


r 1 

2 Vx 

0 


for 0 < x < 1 
else 


Add the graph of the density function using the Nspire-menus given below. The function can be edited 
by using the respective template for a piecewise defined function. 



Maths Templates 

,0.4685 

vise function... 

1 _ . „ . „ 

Double-click on icon to insert item 

% a Vn Vd e D up {" 
{l |o| ooo [S8] M [□] [III] 

loo 

□ 

S n 

a=n- 

t 

piece piecev 



You are invited to double-check in the notes application that the (improper and defined) integral of the 
density function over the interval [0, 1] gives the expected result 1. (You can do this also in a 
calculator page.) 

’l 



dx ► 1 


0 



D-N-L#88 


Guido Herweyers: Statistics with Tl-Nspire 3.2 (4) 


P 7 


(3) Functions of two random variables 


Take any random number X and another one 7 both from the interval [0, 1]. Assume that both are uni- 
formly distributed in [0, 1]. 

Using these two variables new random variables can be defined, as for example: 

S = X+Y,P = X Y,M= max(x, y\ K = min(X, 7), etc. 

What is the probability distribution of these newly generated variables? 

Let’s experiment with random variable S = X + 7 first. Take two random numbers between 0 and 1 
and form the concrete sum s = x+y. Then repeat this experiment very often. 

Numbers x and y are created following the uniform distribution over [1, 0]. Then the sums are values 
between 0 and 2. The density diagram of the generated sums shall become closer and closer to the 
density function of S. 

Do the same for variables P, M and K. 

sample_size n:=3 ► 3 

x:=rand(n):y:=rand(n) (Two commands are separated by a colon. They are both executed when pressing 
ENTER.) 

x ► {0.280027,0.588786,0.734613} 
y ► {0.166901,0.023942,0.38266} 
summ:=x+y ► { 0.446928, 0.612728, 1.11727} 
prod:=x* y ► {0.046737,0.014097,0.281107} 
maxi:=max(x,y) ► { 0.280027,0.588786,0.734613} 
mini:=min(x,y) ► {o. 166901,0.023942,0.38266} | 

f" B 

This is the density histogram for a sample size n = 2500 together with the respective probability 
density function. 



[1] Information how to find the probability density function f4(x) (pdf) is given later in this DNL. 
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Exercise: 

Take random variable K = min(X, Y) from the last example. 

Simulate a sample of size 2500 and find the sample mean x. This is an estimation of the theoretical 

i 

expected value E(K ) = ju = Jx • (2 - 2x) dx. 

o 

Find also the sample standard deviation (stdevsamp), which is an estimation for the theoretical stan- 


dard deviation a - ^/Var(AT) = j(x- jd) 1 • (2 - 2x) dx. 


Calculate fi and a. Check if sample mean and sample standard deviation in fact are lying close to // 
and a. Repeat the experiment for some samples of size 2500 or larger. Repeat the experiment for some 
samples of size 10. What is your conclusion? 


Part 4: 

For the statistics background see [6], [7], [8], [9], [10], [11]. 

Example 1: Does the student say the truth? 

A test consists of 10 multiple choice questions with four possible answers one of them being the cor- 
rect one. A correct answer gives one point, a false one gains no point. 

A student reaches 6 points of the 10 possible ones and he assures that he had made a random choice 
for each question (because he didn’t study the subject). Can you believe him? 

By simulating this test and keeping the scores by automatically data storing (capture) one can investi- 
gate how often a score of 6 or higher is appearing. 


answcr:~{ "correct", "false'V'false", "false"} * { "correct", "false", "false", "false"} 
test * { "false " , "correct " , "false" , "false" , "false " , "false " , "false" , "correct" , "false " , "false " } 

The result of this test is: 
result: -countlfitcst/ - "correct") ► 2 
number; =diml score) ► 2055 
frequ -countIf(score ? ?£6) ► 47 

1.* f requ 

The relative frequency (by experiment) for getting 6 or more points is ► 0.022871 

number 


|®test 

■ score 

■ 

♦ =randsamp(answer,10) 

=capture(result,1) 


1 false 

3 


2 correct 

2 


3 false 

0 


| false 

3 


5 false 

1 


6 false 

■ 

2 


7 false 

5 


8 correct 

3 


B . . 



A4 =" false" 




I 0 


3 4 5 6 

score 


8 9 10 


p 10 
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The experimental frequency of gaining a score above 5 is according to a simulation of 2055 tests only 
-2.3%. It seems to be sure that the student does not tell the truth. 

Let X the correct answers when guessing at each problem. Then X follows a binomial distribution with 
yi — 10 and success probability p = 0.25. 

Let’s check the experimental frequency table by comparing it with the theoretical probability distribu- 
tion: 



The hypothesis test is performed as follows: 

The null hypothesis is set: H 0 : p = 0.25 (the student guessed) 

The alternative hypothesis is: Hi ■ P > o .25 (the student studied and did not guess with each 

question having the same success probability) 

The test variable is the number of correct answers X with X ~ 5(10,0.25) - assuming that H 0 is true. 

The observed value of the test variable is x = 6. The exceeding probability or the p-value is 
P(X> 6) = P(6 <X< 10) = binomcdf(1 0,0.25,6, 10) = 1.97%. 


binomCdf 

i ' 

10,— ,6,10 
4 , 

j 0.019728 


l-binomCdf| 

K 5 ) 

0.019728 


This probability is very small, thus H 0 is rejected which means that we can be (almost) sure that the 
student does not tell the truth (and he did study). 
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Example 2: Is this die a correct one? 

When rolling a correct die 120 times one can expect that the numbers 1 through 6 are appearing ap- 
proximately 20 times each. Which deviations of this expectation are possible? A measure for the de- 
viation of the expected values is the 


Chi-Square value: j 2 = ^ 


(observed value - expected value) 2 
expected value 


=Z 




If this random variable turns out to be small then this die should be correct. But if one finds a large 
X 2 -value then it is possible that the true probability distribution does not correspond with the expected 
discrete uniform distribution. 

For getting an impression of the distribution of the random variable % 2 it will be sufficient to perform a 
large number of simulations in order to create a density diagram of the observed % 2 -values which will 
become stable on the long run. 


The distribution is right skewed. As you can learn from the graphic representation the % 2 -distribution 
with 5 degrees of freedom forms a good model for the probability distribution. 



Roll a die 120 times in order to check if it is a correct one. Assume that you will get the following 
results: 


numbers 

i 

2 

3 

4 

5 

6 

observed 

17 

12 

23 

18 

25 

15 

expected 

20 

20 

20 

20 

20 

20 
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The respective % 2 -value can be calculated as: 

2 (17 -20) 2 (12-20) 2 (23-20) 2 (18-20) 2 (25-20) 2 (15-20) 2 

Z ~ 20 20 20 20 20 20 

Can we expect this result? The /rvalue is (we will ask TI-NspireCAS) . . . 


■ 6.8 


Calculator 

* 

^Actions 

► 

^►.5 2: Number 

► 

X= 3:Algebra 
fd 4: Calculus 

► 

5: Probability 

► 

X §: Statistics 

► 

[□□] 7: Matrix & Vector 

► 

$ € 8: Finance 

► 

ioi ^Functions & Programs 



1:Stat Calculations 
2:Stat Results 
3: List Maths 
4: List Operations 


6: Confidence Internals 
7:StatTests 


1:Normal Pdf... 
2:Normal Cdf. .. 
3:lnverse Normal.. 
4:t Pdf... 

5:t Cdf... 

6:lnverse t. .. 

7:x J Pdf. . . 

9: Inverse x 2 -- 



X 3 Cdf(6.8 ? <»,5l ► 0.235945 


. . . 0.236. Thus, the chance that the value of the % 2 -statistics is at least as great as the observed value of 
6.8 is 23.6% and this is not extra ordinary. 

The null hypothesis can not be rejected by this “goodness-of-fit test”: 

H 0 : the die follows a uniform discrete distribution 
versus 

Hi: the rolled numbers are not uniformly distributed. 


The question is: which observed % 2 -value will lead us to a rejection of the null hypothesis at a signifi- 
cance level of a = 5%? 

This critical value is obtained by using the inverse of the % 2 distribution function. 



The critical value is 1 1.07. 

The critical value can also be found using a slider: 
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Example 3: Testing a proportion 

A sweets manufacturer affirms that 30% of a special sort of candies is yellow. You want to check this. 
So you buy a bag of 50 sweets and find out that 23 of them are yellow - this is 46% ! 

Is this fact sufficient for concluding that the manufacturer does not tell the truth? 

Let p the proportion of the population. 

The hypothesis test is H 0 : p = 0.3 versus Hi \p± 0.3 (a two-tailed test). 


The simulation: 

Assume that H 0 is true. Take a sample of 50 candies out from a population containing 30% yellow 
sweets. 

Investigate how much the sample proportion p can diverge from the (fixed) population proportion p. 
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H 1 


sample 

| ► { "other", "other", "other", "yellow", "other", "other", "yellow", "other", "other", "other", "other", "other", "other", "other", "yellow' 

counties ampl e, ? = " y e 1 low ' ' ) 

p_rand:= ► 0.36 

n 

samp_numb : =dim(samp_rn d) - 1501 ' 

fr e q : =countI f( s amp_r nd, ? > 0. 46 ) ► 23 


experimental two-tailed rejection bound: 


2 freq 
samp_numb 


0.030646 



A sample 

B samp_rnd c 

♦ 

=randsam| 

=capture(p_ranc 

1 

other 

0.28 

2 

other 

0.34 

3 

other 

0.38 

4 

yellow 

0.2 

5 

other 

0.34 

6 

other 

0.24 

7 

yellow 

0.46 

8 

other 

0.34 

9 

other 

0.18 

10 

other 

0.3 

11 

other 

0.26 

12. 


n a 

A7 =" yellow" 


TTT 

6 . 0 - 

5.0- 

4.0- 

3.0- 

2 . 0 - 
1 . 0 - 


0 . 0 - 


0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 

samp_rnd 


The experiment shows that the two tailed exceeding probability is 3.06% (having drawn 1501 
samples). The conclusion is that we can reject H 0 with a 5% significance level. 

Although the probability distribution of the sample proportions is a discrete one - which will become 
clear taking a class width less 0.02 - it can be approximated by a normal distribution with mean 

/ 0.3-0.7 

u = 0.3 and standard deviation cr = J « 0.65. 

V 50 
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Taking the normal distribution as probability model the two tailed exceeding probability is given by 
2 • P(P > 0.46) « 1 .4% (without continuity correction!) 


0.3- 0.7 


0.064807 


50 


2* normCdf| 0.46, 1,0.3, 


0.3- 0.7 


50 


0.013555 


The prop z Test with 1 proportion works with the standard normal distributed test variable 

v P-0.3 u ^ • 0.46-0.3 

Z - , =. The observed value is z = , « 2.47. 


0.3 -0.7 


0.3 -0.7 


50 V 50 

Then the exceeding probability or the p-v alue is 2 P(Z > 2.47) ~ 1.4%. 
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Comment of the Editor: 

Translating and working through Guido’s Cahier was really a pleasure. I learned a lot about 

statistics and a lot about Tl-Nspire’s latest version. 

Two open questions remained: 

(1 ) How to find the probability density functions for functions of random variables 
(Paragraph 3)? In my books I found the convolution integral for the sum X + Y but 
that was all. I asked Guido and he promptly sent the required information, thanks for 
this. You can inform in an extra contribution in this DNL. 

(2) I wondered if it were possible to perform Guido’s experiments and simulations with 
DERIVE, too with not too much efforts. As you will see on the following pages, it 
worked. I used the opportunity to write a few short functions/programs to implement 
some Nspire-functionalties for DERIVE like randSamp in DNL#87 (in most cases 
earlier I did the other way round: DERIVE tools for Tl-Nspire). The functions for 
plotting the diagrams have been developed in an earlier Statistics-Tool-contribution. 
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Following Guido’s Cahier I start comparing the distribution of z-scores and t-scores. 

I implement randnormCmean, standard deviation, sample size) for DERIVE. The counting 
loop is necessary for avoiding repeated samples. 

randnorm(x_, s_, ss, dummy, i) := 

Prog 

dummy := RANDOM (0) 

i := 0 

#1 : Loop 

i : + 1 

If i = 1000 exit 

VECTOR (RAND0M_N0RMAL(s_, x_) , k, ss) 

#2: samples(x_, s_, ss, n) := VECTOR (randnorm(x_, s_, ss) , k, n) 

#3: samplesClOO, 10, 4, 10) 

"111.1 113.9 S3 . 86 130.3 

83.49 98. 68 92.46 S2.44 

96.36 97. SS 10S.4 91.39 

99. 87 102. S 94.33 103.3 

87 . 34 94.62 102 94.73 

#4: 

103.3 S3 . 1 104.6 109 

121.2 101.6 78. 23 107.1 

107. S 93.31 114.4 113.1 

99.32 93.44 107.6 93.06 

_ 109.4 116.6 91.84 103.4 

z_t_scores(mean , stdev, sample size, number of simulations) produces two lists of 
corresponding scores (according to the columns presented on page 3). 

z_t_scores(x_, s_, ss, n, j, sd, smp, sm) := 

Prog 
J := 1 

zscores := [] 
tscores := [] 

Loop 

If j > n exit 

#3: smp := APPRGX(randnorm(x_, s_, ss)) 

sm := AVERAGE (smp) 
sd := STDEV(smp) 

tscores := APPEND(tscores , [(sm - x_)/(sd/-/ss) ]) 
zscores := APPEND(zscores , [(sm - x_)/(s_/-/ss) ]) 

J :+ 1 

RETURN M zscores and tscores stored 11 

The scores are stored in two lists generated of 1500 simulations. Comparing max and min 
values shows the different dispersion. We cannot perform it dynamically like with Nspire. 


#6: 

z_t_scores(100 

, 10, 4, 13O0) 



#7: 

zscores and ts< 

cores stored 



#8: 

[MIN(zscores) , 

MAX(zscores)] = 

[-2.95, 

3.814] 

#9: 

[MIN(tscores) , 

MAX(tscores) ] = 

[-16.09, 

15.25] 
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Like in most cases the graphs are much more impressive. (normpdf(x,g,c>) is not implem- 
ented in DERIVE, normpdf(x) is the N(0,1) distribution). 


#12: 

HI STDENS (zscores , 

-4, 4, 16) 

#13: 

normpdf (x , p := 0 , 

cr := 1) := - 

1 

#14: 

TABLE(normpdf 0O , 

x, -5, 5, 


1 

. F V P 

r 1 

^ x - p ^ 

2^ 

tr.VC2.7r) 

^ 2 

^ cr j 

j 



The TABLE-function is used to obtain a thick Gaussian bell shaped curve. For finding the pdf 
of the Student t-distribution I define the respective function tpdf(x, degrees of freedom). In 
DERIVE the cdf (the distribution function is implemented as STUDENT(x,n).) 


#15: HISTDENS(tscores, -17, 17, 68) 


#16: tpdf (x , n) := 


7(n-rr) -f 


#17: TABLE (tpdf (x, 3), x, -20, 20, 0.01) 



#18: z_t_scores(50 , 2, 10, 1500) 

#19: zscores and tscores stored 

#20: [MIN(zscores) , MAX(zscores)] = [-3.119, 2.901] 
#21: [MIN(tscores) , MAX(tscores) ] = [-3.747, 4.261] 
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#22: HISTDENSCzscores, -4, 4, 32) 

#23: TABLE (no rmpdf 00 , x, -5, 5, 0.01) 



#24: HISTDENS(tscores, -4, 5, 36) 

#25: TABLE (tpdfO, 9), x, -5, 5, 0.01) 



I proceed with the pdfs of one random variable, rand(n) generates a sequence of n uniformly 
distributed random numbers from [0,1] (with n > 1). For n = 1 take RANDOM (1). 

#26: rand (n) := VECTOR (RANDOM CD, k, n) 

#27: rand(5) = [0.67S7, 0.S265, 0.4101, 0.5258, 0.4188] 

#28: HISTDENS(rand(200), 0, 1, 10) 



I start with sample size 200 and increase up to a sample size of 10 000 which wouldn’t be so 
easy done with Nspire! 
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HI STDENS( rand (1000), 0, 1, IQ) 



HISTDENS(rand(4000), 0, 1, 10) 



HISTDENSCrandCLOOOO), 0, 1, 10) 



Now we generate 200 random variables x and study the distribution of the squares. 


2 


HI STDENS (VECTOR ( rand (1) , k, 200) , 0, 1, 10) 



HI STDENS (VECTOR ( rand (1) , k, 2000) , 0, 1, 10) 



Can you guess the form of the probability density function? 
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2 

#34: HI STDENS (VECTOR ( rand (1) , k, 5000), 0, 1, 10) 


#35: 


1 


TABLE 




x, 0, 2, 0.001 



What about U = 3X 3 ? 



I superimposed two 5000 simulations runs and then added the respective pdf. I am quite 
sure that you can easily follow how to derive the density function (compare with page 6). 

The density function of U is the first derivative of the distribution function TV of U: 
For0<x< 1: 


f- 

h 

F u (x) = P{U <x) = P(3X 3 <x) = P(X < 3 -) = J 1 • dt 


1 = 3 — 


Hence, the density function of U is given by: f v (x) = FJ (x) = < 


1 

V 81 x" 

0 


for 0 < x < 1 

. Does it fit? 


else 
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Functions of two Random Variables 

Distribution ofX+ 7 with 3000 simulations performed: 


#39: HISTDENS (VECTOR (rand Cl) + rand(l), k, 3000), 0, 2, 40) 

fCx) := 

If x < 0 
0 

If x < 1 
#40: x 

If x < 2 
-x + 2 
0 

#41: TABLE(f (x) , x, 0, 2, 0.001) 



Compare with the pdf of U= 2X. 


#42: HISTDENS (VECTOR (2 ■ rand(l) , k, 3000), 0, 2, 40) 

f 1 1 

#43: TABLE , x, 0, 2, 0.001 

l 2 J 



It should be no problem to find the pdf flu) = 0.5. 
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Part 4 


Example 1: 


I need my DERIVE function randsampCpopul 

randsamp(l, n, s := 0, dummy, i, samp) : 
Prog 

dummy := RANDOM (0) 

i := 1 

samp := [] 

Loop 

#44: If i > n 

RETURN samp 

k := RAND0M(DIM(1)) + 1 
samp := APPENDfsamp, [lj,k]) 

If s = 1 

1 := DELETEfl, k) 

i :+ 1 

#45: ans := [1, 0, 0, 0] 


#47: results(ans, 10, 2000) = results 


ion list, sample si ze), s = 0 by default. 

results(l ist , k, n, t, i_, j_) := 

Prog 

i_ := 0 

res := [] 

Loop 

If i_ > n 

RETURN "results stored in res 11 
#46: t := I(randsamp(l ist , k)) 

res := APPEND(res, [t]) 

J- := 0 
Loop 
J- : + 1 

If j_ = 500 exit 

i_ :+ 1 


red in res 


This gives the number of correct answers simulating 2000 tests consisting of 10 multiple 
choice questions each. 


FRETAB(list) and FREQDIAGO i st) as well are functions provided in my statistics tools utility 
file. If you load the stat_4 . dfw file all functions and programs needed are available and 
you can see them via Author > Function Definition. 


#48: 


#49: 


# 50 : 


FRETABCres) 

'01 2 3 4 5 67' 

. 88 389 574 511 288 112 27 12 _ 

FREQDIAG(res) 



We have 88 tests with no single correct answer and there at the other side 12 tests with 7 
correct answers. There are no tests with 8 or more correct answers - with randomly chosen 
answers!! 

HISTO(list, start, end, number of equal wide cl asses) plots the respective histogram, 
#51: HIST0(res, -0.5, 10.5, 11) 
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Example finished: 


#52: 

#53: 

#54: 

#55: 

#56: 


DIM(5ELECT(k > 6, k, res)) = 39 


39 

= 0.0195 

2000 


r n 

1 - BINOMIAL_DISTRIBUTION 5 , 10 p 

l 4 


0.01972 


res ults Cans , 10, 5000) = results stored in res 


DIMCSELECTCk > 6, k, res)) 
5000 


0.0188 


Example 2: 

randint(a,b,n) generates a? integer random numbers x with a<x<b. si mul (number) gen- 
erates number die rolls. 

#57: randint(a_, b_, n_) := VECTOR (RANDOM(b_ - a_ + 1) + 1, k, n_) 

#58: randintflp 6, 20) 

#59: s imul (n_) := randint(l, 6, n_) 

Make a first test with 120 experiments: 


#60: 

test := s inm"! (120) 


#61: 

test := [2. 3, 3, 2, 2, 3, 3, 2, 

2, 3, 2, 

#62: 

FRETAB(test) = 

'12 3 4 

. 9 16 27 21 

5 6 ' 

25 22 . 

#63: 

obs_ := (FRETAB(test)) 

2 


#64: 

exp_ := [20, 20, 

20, 20, 20, 20] 



I start performing the x 2 -test (same variable names are used as in Guido’s paper): 

2 

(obs - exp ) 

DIM(obs) i i 

#65: chisqu(obs, exp) := I 

i=l exp 


chisqus (n , i , s im, o) := 

Prog 
i := 0 

chiqs := [] 

Loop 

#66: If i = n 

RETURN Ir cbisquares stored in chiqs" 
s im := s imul (120) 

0 :z VECTOR(FREQ(a, sim), a, 6) 

chiqs := APPEND(chiqs , [chisqu(o, exp_)]) 

1 :+ 1 


#67: chisqus(2QQ0) 
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#68: chisquares stored in chiqs 

#69: [MIN(chiqs), MAX(chiqs)] = [0.1, 21.3] 

#70: HISTDENSCchiqs, 0, 22, 22) 


Before plotting the density diagram I inform myself in the DERIVE Online Help if the 
% 2 -pdf is available or not? It is not. What I find is CHi_SQUARE(x , n) which is the x 2 -cdf - the 
distribution function with variable x and n degrees of freedom. 

I define my own x 2 -pdf as follows and proceed according Guido’s guide line. 


1 df/2 - 1 - x/2 

X2_pdf(x, df) := -x -e 

#71: ' df/2 f df 


#72: x2_pdf(x> 5) 

- x/2 3/2 

^2 ■ e -x 
#73: 

6 ■ JW 

#74: TABLE Cx2_pdfCx, 5), x, 0, 22, 0.01) 



The result is a really nice and convincing plot, isn’t it? 

We could do without defining the pdf remembering that the pdf is the derivative of the cdf 
(probability distribution function): 

Please compare #73 from above and #77 on the next page! 
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d 

#75: — CHI_SQUARE(x, 5) 

dx 

- x/2 2 - x/2 

v f 2-e ■ (x + 3-x + 3) ^2-e -Cx + 1) 

#76: 

6 • Jn ■ ,/x 2 ■ JW ■ Jx 


- x/2 3/2 

V2 ■ e -x 

#77: 

6 -Jn 

Calculating the critical value is an easy job now: 

#78: 1 - CHI_SQUARE(6 . 8 , 5) = 0.2358 

#79: N SO L VE ( CH I_SQU AR E ( x , 5) = 0.95, x) = (x = 11.07) 

#80: N50LUTI0NSC1 - CHI_SQUARE(x, 5) = 0.05, x) = [11.07] 


Working with a slider is very attracting with Tl-Nspire. With little fantasy we can do a similar 
animation with DERIVE, too: 

We plot the x 2 -pdf and the (red) segment - expression #81 - with its point ( 1 , 0 ) on the x-axis 
moveable. #82 is the area under the pdf for x > /. 


#81: 


1.5 - 0.5-1 

1 0.1329-1 -e 


#82: 1 - CHI_SQUARE(1, 5) 



The value of the area is given by the distance between the blue line and the x-axis. Moving 
the slider moves the segment which moves the blue line. Try to move the blue line to a 
distance of a = 5% and then read off the value of /. This is it. 


Note: Finding the pdf of the t-distribution as derivative of student (x , n) is possible but not so 
easy. 
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Example 3: 

I believe that the following expressions are self explanatory. 

#83: population := [yellow, yellow, yellow, other, other, other, other, other, other, other] 

#84: sample := rands amp (population , 50) 

DIM(SELECT(c = yellow, c, sample)) 

#85 : p_rand := 

DIM(sample) 

#86: samp_rnd(n) := VECTOR(p_rand , k, n) 

2 ■ DIM(SELECT(v > 0.46, v, samp_rnd(Z000))) 

#37: = 0.023 

2000 

We add the plots: 


#88: HISTDENS(samp_rnd(2Q0Q), 0, 1, 50) 

r r 0.3 -0.7 y 

#89: normpdf^x, 0.3, ^ J 



And this is the remaining calculation done with DERIVE: 


r r ro.3-o.7y r r 0.3-0.; 

■ iNORMALll, 0.3, A - NORMAL lo . 46, 0.3, A 

r r ro. 3 -o. 7 y r ro.3-c 

-INORMALll, 0.3, A II - NORMAL jo . 46, 0.3, A 


f5.j21'| rs.721 1 ! 

ERF l— J - ERF hrJ 


8-^21 


0.01355 


0.46 - 0.3 


•ra 


2.468 


2- (1 - NORMAL (2. 468)) = 0.01358 
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Probability Density Functions (pdf) for Combined Random Variables 

Guido Herweyers, Belgium 


Let X and Y be independent random variables with a common uniform distribution on the interval 
[ 0 , 1 ], having p.d.f. f (x) = 1 for 0 < x < 1 and 0 otherwise. 

Then the joint p.d.f. of X and Y is g(x, y) = f(x)f(y) = 1 for 0 < x < 1 , 0 < y < 1 and 0 otherwise. 


1) Distribution of S = X + Y 

Let S = X+ Y, then the c.d.f. of S is H(s) = P(S < s) = P(X + 7 < s) = II g(x,y^dxdy , 

A 

where A = {(x, j)|x + 7 < s } . 

0 , 5 < 0 

Let R be the square R = [0, l] x [0, l] , then H(s ) = < jj dxdy = area (A nR) ,0<s <2 . 

AnR 

1 , s>2 




0 , s < 0 
5 , 0< s <1 
2 - 5,1 < s < 2 
0 , s >2 


The p.d.f. of S is /z (^v) = 


dH(s) 


ds 


Remark: 

The distribution of S = X + Y , for independent random variables X and 7, is the convolution of the 
distributions of X and 7: 

oo 

(f*gXs) = h(s)= j f(s-x)-g(x)dx 


f and g are the densities of X and 7 respectively. 
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This is done as follows: 

E/(0,1), Y ~ E/(0,1) <mdS = X+ Y with 0 < w < 2 

0 ,s < 0 

w w 

| / MgC? _ x )dx = ^dx = w, 0 < w < 1 


h { s ) = 


i i 

| / ( x )s( s _ x )dx = | dx = 2 - w, l<w<2 


W — 1 

0 ,s >2 


2) Distribution of S = X • Y 


Let 5 = X • Y, then the c.d.f. of S is H(s) = P (S < s) = P (X • 7 < s) = JJ g (x,y)dxdy , 

A 

where A = {(x, y) |x • y < 5 } . 


Let R be the square R = [0, l] x [0, l] , then H(s ) = 


0, 5 < 0 

|| dx dy = area (A nR) , 0 < s <1 


1, s > 1 



The p.d.f. of S is /z(s) 


dH(s) 

ds 


0, s < 0 
ln(s), 0 < s < 1 
0, s>\ 


The value h{ 0)can be assigned arbitrarily, e.g h( 0) = 0. 
Remark: the p.d.f of S is not bounded on the interval [0,1] ! 
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3) Distribution of M — max(X,r) 

Let M= ma x(X,7), then the c.d.f. of Mis H{m ) = P(M < m) = P(max(X,7) < m) = || g(x,y)dxdy , 

A 

where A = {(x, y) | max(x, y) < m} = {(x, y) \ x < m and y < mj . 


Let R be the square R = [0,1] x [0,1] , then H{m) = 


0, in < 0 

|| dx dy = area (^4 ni?),0<m<l, 


1 , m > 1 



area H (/??) = /?7 2 ,0 < m < 1 


The p.d.f. of Mis /z(m) 


dH[m) 

dm 


0, m < 0 

2m, 0 < m < 1 

0 , p> 1 


The value /?(1) can be assigned arbitrarily, e.g. h( 1) = 2. 


4) Distribution of K = min(A",T) 


Let K = min (X,Y), then the c.d.f. of K is H(k) = P (K <k) = P(min(X, Y)<k) = JJ g(x,y)dxdy, 

A 

where A = {(x,y) | min(x,y) < k} = {(x,y) \ x< k and y^k}. 


Let R be the square R = [0,1] x [0,1] , then H(k ) = 


0, k< 0 

JJ dx dy = area(ylni?) , 0 < k < 1 , 


k> 1 


The p.d.f. of K is 


h(k ) 


dH(k ) 

dk 


0 , k < 0 

< 2-2 k, 0<k<\ 
0, k > 1 


y 



The value h( 0) can be assigned arbitrarily, 
e.g. h( 0) = 2. 
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I must admit that this was completely new for me and I - nasty me - sent more problems to 
Guido and asked for some reference literature because the only facts which I could find were 
the notes about applying convolution for X + Y. They are given on page 29. 

I performed some simulations using TI-NspireCAS and DERIVE as well and wanted to find 
the formulae for the respective density functions. My problems were: 

5) Distribution of S = 5X -2Y 


HISTDENS(VECT0R(5 • RANDOM(l) - 2.RAND0M(1), k, 5000), -3, 6, 45) 




Inspecting the density diagrams I had the impression of a trapezoidal distribution and I tried 
to find the boundary lines as you can see left above. 

And this is how Guido treated this distribution: 


Distribution of S = 5X - 2 7 


The reader can verify that the p.d.f. of S is h (s) 


dH(s) 

ds 


0, s < -2 
0. Is + 0.2, — 2 < s < 0 

< 0.2, 0<s<3 

— 0. Is + 0.5, 3<s<5 

0, s > 5 


Remark: Let S = X — Y , Xhas a uniform distribution on the interval [0,5] and 7 a uniform distribu- 
tion on the interval [0,2], then we expect that S has the same distribution as problem (6). 

Indeed, the joint p.d.f. of X and 7 is g (x, y ) = — • — = 0. 1 for 0 < x < 5 , 0 < y < 2 and 0 otherwise. 

The p.d.f. of S is (see M. H.Degroot, M. J. Schervish, Probability and Statistics , Fourth Edition, 
Pearson International Edition, 2010, page 178): 


0, s < -2 

2 

1 0. 1 dy - 0.2 + 0. Is, - 2 < s < 0 


00 

H 5 )= 1 s(s + y,y)dy- 


A 

Jo.lJy = 0.2, 0<s<3 


| O.ldy = 0.5-0. Is, 3<s<5 

0 

0, s > 5 


This is exactly the trapezium from above! My conjecture is confirmed. 
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For my next problem I was not able to find the density function by inspection only. 


6) Distribution of S = X 1 2 + Y 2 

2 2 

HISTDENS (VECTOR (rand (1) + rand(l) , k, 5000), 0, 2, 40) 



All what I could imagine was that the first part (0 < x < 1 ) of the density function seems to be 
a horizontal line (y « 0.8?), and what’s the second part?? 

Here is Guido’s answer: 


Distribution of S = X 2 + Y 2 


Let S = X 2 + Y 2 , then the c.d.f. of S is H(s) = P(S < s) = p(x 2 + Y 2 < 5 ) = jj g{x,y)dxdy , 

A 

where A = j(x, y ) | x 2 + y 2 < . 

Let R be the square R = [ 0 , l] X [ 0 , l] , then 


H(s) = 


0 , 5 < 0 

|| dxdy = area (A nR) , 0 < s < 2 = 

AnR 

1 , s>2 


0, s < 0 
ns 

— , 0 < s < 1 

4 


1- | s - x 2 ^jdx ,1 < s <2 

A7A 


1 , 5 >2 


The p.d.f. of S is h(s) = 


dH(s ) 


ds 


0, s <0 


K 


0 < 5 < 1 


n 

arctan 

4 


(V^T),1<5<2 
0, s >2 


The value h (0) can be assigned arbitrarily, e.g. h (0) = 0 . 


D 
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DERIVE gives another output for 
the integral. But plotting both 
functions shows the identity for 
both expressions. 

It needs defining the domain for s 
then DERIVE confirms Guido’s 
result. 


H(s ) 


= i-J 


^Cs - 1) 


Cl - Jis - x )) dx 


s ■ATANC-/C S - 1)) 
H(s) := + 


s ■ ASIN 


li) 


+ 7(s - i) 


Ht) 


— H(s) 
ds 


ATAN C/(s - 1)) 


The plot below is really convincing 
(density diagram based on 10 000 
values X 2 + Y 2 ). 


s :e Real (1, 2] 


ASIN 



2 


ATAN C/(s - 1)) Jr 

= ATAN C/(s - 1)) 

2 4 



Guido’ comment: It is interesting to mention that follows a % 2 -distribution with two degrees of 
freedom, if bothXand Y are independent standard normal distributed. (Show this by simulation! !) 

The last example is left for the reader. 

Let’s compare the density functions of V = 5X and U = X+X+X + X + X\ 

HISTDENS (VECTOR (rand (1) + rand(l) + rand(l) + rand(l) + rand(l), k, 5000), -3, 6, 45) 
HISTDENS (VECTOR (5. rand (1), k, 5000), -3, 6, 45) 



Which one is the blue one, and which one is the red one ( U or V)? 
What makes the difference? 

What are the density functions? 
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Von: Robert SETIF [robert.setif@gmail.com] 

An: Josef Bohm 

Betreff: variance of list with weights in Derive 
Dear Josef, 

It seems to me that there is no command for variance of lists of numbers with their weights. 
For instance [61,64,67,70,73] with [5,18,42,27,8] which is 5 times 61,18 times 64,... 

Matlab and XCAS find 8.5275, but Maple_16 finds 12.0207. 

??? 

No command in Mathematica_6, nor in MuPad, nor in Scilab and nor perhaps in XMaxima. 
Best regards. 

Dear Robert, 

Please have a look, this is my DERVE work sheet dealing with your question: 

#1: [numbs := [61, 64, 67, 70, 73], freq := [5, 18, 42, 27, 8]] 


#3: vari (numbs, freq) = 8.5275 

#4: This meets the XCAS/Mat lab-result. 

I repeat with a shorter list given in two ways: 

#5: [til := [2, 3, 5], fl := [3, 1, 4]] 

#6: list := [2, 2, 2, 3, 5, 5, 5, 5] 

#7: vari Cnl, fl) = 1.984375 

DERIVE's variance and standard deviation are the statistics (sample) values = n-1 weighted, 
vari 1 is another form of vari 

2 

#9: STDEV(list) = 2.267857142 


DIMCn) 2 
I n -f 


#2: 



#10: 


DIMCn) 

I 

i=l 



varilCn, f, av) := 


Z(f) 


#11: varilCnumbs , freq) = 8.5275 

vari 2 is the DERIVE way. 


#13: 



vari2Cn, f) := 


ZCf) - 1 


#14: vari2Cnl, fl) = 2.267857142 

#15: VARIANCES is t) = 2.267857142 

#16: vari2Cnumbs , freq) = 8.613636363 
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But this does also not look like the Maple result. 
And this is the third way to achieve the same result: 


1 DIM(rt) 2 

#19: vari3(n, f) := ■ I n -f 

XCf) i=l i i 



2 


#20: variBCnl, fl) = 1.984375 

#21: variBf numbs, freq) = 8.5275 


So you can be quite sure that 8.5275 is the correct result. I have no idea where the 
Maple-result comes from. 

To be on the safe side, I repeated the calculation with the Tl-Nspire statistics tool 
which confirms the earlier results. 



Best regards 
Josef 
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Giuseppe Ornaghi [g.ornaghi2@tin.it] 

Is it possible, and how, to define in Derive LI function (logarithmic integral) for a complex 
number? 

Thank you very much, 

Giuseppe 


Fred Tydeman’s problem - A non recursive definition of a special sequence 

In DNL#87 Fred complained that he only received answers from MATHEMATICA and Maple 
Users. Fortunately DUG members also tackled the problem: 

This is Benno Grabinger’s answer: 

The number of fix point free permutations of n elements is given by 

a(l) = 0 

a(n) = n ■ a(n - 1) + (-1)" 

It is easy to derive the explicit formula from the recursive definition: 
a(l) = 0 

a(2) = 2a(l) + l |:2! 

a(3) = 3a(2)-l |:3! 

a(4) = 4a(3) + l |: 4! 

a(n) = na{n-\) + {-\) n | :n\ 

This gives: 

a(l) = 0 

a (2) _ q(l) J_ 

2! “ 1! 2! 
a (3) a (2) 1 

3! “ 2! 3! 

a( 4) _ a(3) 1 


a(n) _ a(n - 1) (-1)" 

n\ («-l)! n\ 

a(ri) _ 1 1 1 _ (-1)” 

n\ ~ 2! 3! + 4! + " n\ 


a(n) = n\ 1 — 

2! 3! 4! 


” ! y 


adding all equations leads to: 


Hence, the desired explicit formula is 


a(n) = n 

k=0 


(-1) A 

k\ 
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Substituting in the power series e x = V — for x = -1 then 



On a problem by Fred Tydeman 


Stefan Welke, stefanwelke@web.de 


In DNL #87 Fred Tydeman asked for a non recursive definition of the recursively defined 
sequence / (0) := 1 and f(n):=n • f(n-\) + (-1)" for n > 0 , which is closely related to the 
factorial function, see the remark at the end. The computation should be done with DERIVE. 
This definition is easily seen to be equivalent to the following two step recursion, which has 
some advantages: 


Note, that we start here with n = 0 in contrast to Fred. Now we turn this definition into an 
iteration by the observation that the two step recursion above ist equivalent to the following 
matrix equation: 


/ (0) := 1, / (1) := 0, and/O) :=(«-!) (f(n-l) + f(n -2)) forn>l 


( 0 . 1 ) 



n~ 1 n-\)yf(n-\)) y f(n ) J 


We set M k := 


(o n 

and get by iteration: 


( 0 . 2 ) 



n 


Here the product is meant as P| M k = M n ■ M n , • ... • M x , because the matrices M k do not commute 


for different values of k . The equation (1.2) gives a non recursive definition for consecutive pairs of 
elements of Fred’s sequence for arbitrary initial values /( 0) and /( 1) . 


In Fred’s case we have /(!) = !• /(0) - 1 = /(0) - 1 , so we arrive at 



(0.3) 


for n > 0 . 
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A straightforward implementation in DERIVE is the function g , which works with matrix 
multiplication, the initial value for n = 0 is p : 

g(n , p) := 

If n = 0 
P 

Prog 
k := 1 

#1: m := [0, 1; 1, 1] 

Loop 

k := k + 1 
If k > n 

RETURN (m.[p; p - l])j,lil 
m := [0 , 1; k, k] *m 

#2: VECTORCgCq , 1), q, 0, 7) 

#3: [1, 0, 1, 2, 9, 44, 265, 1854] 

#4: VECTORCgCq, 0), q, 0, 10) 

#5: [0, -1, -1, -4, -15, -76, -455, -3186, -25487, -229384, -2293839] 


An even better and faster implementation looks at two consecutive vectors and avoids actu- 
ally the matrix multiplications: 

(7("-l)l ( An) U An) 

y fin) J (f(n + 1)J [n-(f(n-l) + f(n)) 

This reads as a DERIVE function: 

f Cn , p) := 

If n = 0 
P 

Prog 
k := 0 

#6: v := [p, p - 1] 

Loop 

k := k + 1 
If k > n - 1 
RETURN vj,2 

v := [vj.2, k-Cvj.1 + vj.2)] 

#7: VECTORCfCr, 1), r, 0, 15) 

#8: [1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961, 14684570, 

176214841, 2290792932, 32071101049, 481066515734] 

And in 0.030 seconds on my quad-core computer: 


x 

for n > 0 

J 


#9: f C1000 , 3) 
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#10: 

952804520525856627451321481964974025676536975193112726642291666949765841919494112 

3654074043649936771469190498393074595125199264199274530745748888275624174775905326048951 

2818605825668328628235028158861129380513898824291414271642752061375176455335859217179843 

0047917712959117223849859698300434989843430835148458343532497042559094724612062055140061 

0978315986801768791096183442634217492198922981208732028234731570475335628675073790994250 

8875727649114731243454105244908693902862698057664474017776295243350980101824226317390346 

3572002698088459578448187740801226069477070879140617548140726852460925203915519984473771 

9399545837132580960644563545028399963331023121808491402771841379009745212254143840510300 

6694789923523688814802747424087542627389558146387253075616289154444789420104221571265508 

7310620503990472303537105338553902408633211673278396807021386154622670101589942979864048 

5448836064494677803515601109344904598724771483892921318415306238621955407953171903684274 

8564774292875509855087274087291025196660117293707848372593657093042198138917380394110079 

9271921035231381269018568788603781187966173586134639703398516141805567303108810732584407 

7470063807365854600515371697425242865368084262956008405844626962495146170117473602597621 

0863882405329625430185563257576587431662494598489186336712344414519043646298571584399499 

5585415484768055489426473424303010136097860621068702897445838967472097020663391272597658 

2304058376342313430991772914865965597618464937306241026322365044601712806648180182270156 

7833148345030064644891209701382661038880491585883074341911495897609445361802846034613234 

2910296808936520915655170485317404350873465650015092338112581633790783702873568796288080 

8126034990533482473718607751970880012428181918980694189406702875324728431636272269479825 

4522683334457437048916416453039718615017376271682981758710041801548656827276837251789656 

0822334214967612728071207407986701078577759387340199402815725800534178449672827461257277 

8914007758849393689510446063664661812442326908450100905621441832596197383417383411141235 

0836243372003018439173859692513226260242227288258218281193174256854111903986188219385074 

3695550378632545236274910398571040931773645905998535373429245793593749699723266387296496 

6522096514786017905504973784467529046100136404435690046191512560757962999662268766379946 

7599266677222287792956366453891264317223932685390686011191139084749854518135087503539806 

6868621959973870036473108206470890805125591766035651660263166256071859066523494404932873 

9892430338853873103631687347396881324950657650842869854703810748598526517214826649170192 

27944750044815550686001 

f (10000, 3) needs 0.631 seconds. 






p 40 


User Forum 


D-N-L#88 


An equivalent result is the asymptotic behaviour of / : 

(0.8) / ( n ) ~ + /(0) - 1) • n ! for large n 

Finally we consider the derive function h(n) which computes the values of / according to 
formula (1.6) with a = -1 and /(0) = 1 as Fred’s initially posed sequence. 


#11: h(n) := I 


f k 
(- 1 ) 


•, k, 0, n 


k! 


■n! 


#12: VECT0R(h(q) , q, 0, 7) 

#13: [1, 0, 1, 2, 9, 44, 265, 1854] 


Although the function definition for h looks much simpler than the previous definitions for g 
and f , this implementation has the disadvantage that computation time explodes . h (1000) 
needs 2.82 seconds and h(5000) needs 518.8 seconds. (1.6) is a very pleasant formula in the 
eyes of a mathematician, but it is not suitable for programming purposes. 


And there was a third reaction from Australia sent by our DUG member David Halprin: 

Josef 

In December 2011, I sent you this letter with attachment BOEHM.PDF in the hope 
that you would find space in a DNL. The purpose was to elicit reader feedback as 
well as a challenge to see assorted ingenuities. 

I know you have been very busy. 

BTW I spent many hours of many weeks in a fruitless attempt to solve the Fred 
Tydeman recursive problem with true mathematical reasoning. I read that some 
people came up with a Mathematica or Maple 'solution' but that evades the kernel 
of the problem. 

The first definition explicitly defines the series without recourse to recursive meth- 
odology 

The second definition, albeit recursive is insufficient, since it is incomplete. 

I spent my time trying to find an expression that defines the n th term as efficiently 
as Tydeman's first definition, but with its basis in a recursively defined series, as 
one gets with Fibonacci and Lucas types of series. 

All these types of series can be treated under the one generalisation, called a Gen- 
eral Admixture Series, (G.A.S.) for which I wrote a very detailed paper back in 1989 
covering Tribonacci, Quatracci, Pentacci and much higher, for which all could be 
defined with any of the 5 methods below: - 

1) Each term is the sum of the previous 2, 3, 4, 5, 6, terms etc.. 

2) An equation, which determines the (n+l) th term. 

3) The sum to the (n+l) th term. 

4) The limiting value of the ratio of two successive terms. 
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5) The generating function, 
viz. from 2) above 


(a A n - b A n)/(a - b) for Fibonacci-type series, or (a A n + b A n) for Lucas-type series. 

However, I had to give up, sadly. I hope that someone among the readership will 
achieve a real mathematical solution. 


Herzlichst 
David H 

First of all I’d like to apologize for my not publishing David’s paper from 201 1 . It is following. 
David promised to submit an updated version of his Tribonacci-Quattracci... paper men- 
tioned above. He wrote: 


Josef 

Coincidentally I spent a solid 8 hours today re-attempting the Tydeman series, still 
thinking I had cracked it, but alas and alack, not to be as yet. I was constantly re- 
ferring to my methodology in my paper, attached RECURSIV.PDF. Please put me in 
the queue for publishing it, but NOT this copy of today, since it needs some addi- 
tions, for which I have made some notes within. However please read it and enjoy. 
You will see that my approach is completely different from Stefan Welke and Benno 
Grabinger, whose interesting papers you sent me; thanks very much. 

My paper ends with quite a challenge for the readers. Maybe my hoped-for answer 
is in the too-hard basket for today's mathematics???? I would certainly welcome 
your comments. 

Herzlichst 
David H 


DILEMMA AND/OR PARADOXON 

David Halprin 

I have been attending math, seminars at Melbourne University for years. Most of them are 
under the auspices of MUMS (Melbourne University Mathematical Society", whose members 
are an admixture of students, ex-students and interested parties. 

URL: http://www.ms.unimelb.edu.au/~mums/seminars/pastseminars.html 

These daytime seminars are slotted in between formal lectures and are allowed one hour. 
Recently, there was some time to spare, so the lecturer presented a question to the atten- 
dees for them to come up, one at a time, to chalk their solution(s) on the green board. 
Belatedly, I realised that I should copy down all these solutions, but I only managed to copy 
five; the last two I missed out and could not duplicate, so I am requesting the readers to 
submit their opinions of what the last two solutions were. 

(What were the answers of student 6) and 7)? You will find the solution of this problem in the 
next DNL - provided by David, Josef.) 

(If you are not familiar with D 1/2 which is the 0.5 derivative then read David’s contribution “On the 
lighter Side of Operational Calculus” in DNL#35. Josef) 
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QUESTION 

Solve for y = f(0) 

1 ) First-Year Student 

y = sin0, Z>y = cos0, Z> 2 y = -sin0 Q.E.D. 

2) First-Year Student 

y = cos0, Dy = - sin0, D 2 y = -cos0 Q.E.D. 

3) Second-Year Student 

With respect, you two guys are both correct, but one needs to combine both solutions 
under one function. 

y = c/s0 = cos0 + /-sin0 = e ie 

Dy = i e‘\ D 2 y = -e iQ Q.E.D. 


4) 


Second-Year Student 

With respect, y’all missed out on another avenue of calculus venturing 

Let y = sun0 = Z>^ sin 0, D^sun0 = D^ ■ Z>^sin0 = D sin 0 = cos 0 
Let D ^ cos 0 = cus0, D^cus0 = D 1/1 ■ D' /2 cos0 = D cos 0 = — sin 0 
Z>sun0 = D l/l • Z>^sun0 = D l/l cos0 = cus0 

D 2 sun0 = Z)cus0 = • Z)^cus0 = — D^sin0 = — sun0 Q.E.D. 


5) Second-Year Student 

With great respect, what about 

Let y = cus0 = D l/l cos 0, D^cus0 = D l/l • D l/l cos0 = D cos 0 = - sin 0 
Let D ^ sin 0 = sun0, Z>^sub0 = D'^ • Z)^sin0 = D sin 0 = cos 0 
Dcus0 = D^cus0 = - Z)^sin0 = -sun0 

Z> 2 cus0 = -Z>sun0 = — D l/l ■ Z>^sun0 = —D^c os0 = — cus0 Q.E.D. 

6) Third-Year Student 

With respect, you two guys are both correct, but one needs to combine both solutions 
under one function. 

7) Third-Year Student 

With greatest respect, the combined solution of the third solver and the sixth solver 
need to be combined as one umbral solution, don’t you think? 


David Halprin 
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Titbits (39)- Emulating the EIGamal Cryptosystem (Part 1) 

by Johann Wiesenbauer, Vienna 

Well, there has been quite a break since my last column here due to the fact 
that I was very busy lately. Nevertheless, I'm ready now to continue now with 
one of my favourite topics dealing with the use of elliptic curves in public-key 
cryptography. To be more precise, I would like to make a setup in DERIVE for 
the EIGamal cryptosystem, by far the most important competitor of RSA. Most 
of the routines, which I introduce here, I developed for my talk at DES-TTME- 
conference in Malaga in 2010. If you missed my contribution in the proceedings 
of that conference, this article might make up for it to some extent. 

Let's start with some general remarks on the EIGamal cryptosystem. As you 
might already know, it's based on the discrete logarithm problem or DLP for 
short. In its most general form what you need is a cyclic group G along with a 
generator g such that for any heG the equation g x = h is usually impossible to 
solve in a reasonable time. Needless to say that the group G must be rather big 
to fulfil this condition, otherwise you could find x simply by trial and error. The 
"classical" choice for G is the multiplicative group of the residue class ring 
mod p for some prime p with say at least 1024 bits, but nowadays cyclic groups 
consisting of points on an elliptic curve over some finite field are very popular, 
too. In fact, we are considering here the latter case, where the finite field will 
always be a residue class ring mod p again for some prime p with at least 160 
bits. It turns out then that group G itself has about the same order of 
magnitude, hence it can be much smaller than in the first case, which is one of 
the main reasons of the popularity of E CC (=EI liptic Curve Cryptography) in any 
environment where only a small amount of resources is available, like on 
smartcards or handhelds. 

This said I'll focus more on the mathematics behind EIGamal now rather than on 
other cryptographic aspects. For a start, let me remind you of same facts 
concerning the arithmetic on elliptic curves over a residue class mod p, where p 
is a - usually "big" - prime number (cf. my Titbits 30 if you need to brush up on 
this topic). Basically, we are dealing with an algebraic curve of the form 

y 2 =x 3 + ax + b (a,beZ/pZ) 

along with the point O at infinity, where the discriminant 4a 3 +27b 2 of the 
polynomial on the right-hand side does not vanish mod p, i.e., this polynomial 
hasn't got multiple roots in the residue class ring Z/pZ. Under these conditions 
the finite set E a ,b(Z/pZ) of points on this form an abelian group, if the addition 
of two points U and V on the curve and the n-th additive power of U is defined 
in a way that can be seen by looking at the following DERIVE-routines: 
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add(Up v, a, p, k_) := 

Prog 

If U = [p, p] v v = [p, p] 

RETURN u + v - [p, p] 

If uj,l = v;l a M00(ui2 + vi2, p) = 0 
RETURN [p, p] 

If U : v 

k_ :: MOD((3'Uil A 2 + a).lNVERSE_M0D(2< m2, p), p) 
k_ :: M0D((vi2 - ui2>.INVERSE_HOO(vil - ml, p), p) 
a := M0D(k_ A 2 - uj.1 - v^l, p) 

[a, M0D(k_»(uil - a) - ui2, p)] 

multCu, n, a, p, b_, u_) := 

Prog 

b_ - [p, p] 

u_ := [uil, p - m2] 

If n < 0 

RETURN -n, a, p) 

Loop 

If n = 0 
RETURN b_ 

If M0D(n, 4) = 1 

[n 1, b_ ~ add(b_, u, a, p)] 

If M0D(n, 4) - 3 

[n :+ 1, b tz add(b_, u_, a, p)] 
u := add(u. u r a. p) 
u_ [ml, p - ui2] 
n :/ 2 


Basically, mult(u,n,a,p) makes use of the so-called NAF-representation of n, 
which uses the "digits" -1,0,1 and is unique provided that never two nonzero dig- 
its are adjacent (hence the notation NAF, by the way, meaning "non adjacent 
form"). If n is even then the least significant digit is always 0, otherwise it is 1 
or -1, depending on whether n = 1 mod 4 or n = -1 mod 4. By removing this last 
digit and responding to its value in a similar way as in the original "Square and 
multiply"- method, except for using inverse points whenever -1 occurs, one can 
proceed in an iterative way as can be seen above until the stop condition n = 0 
becomes true. Actually, it is supposed to be a little bit faster than the normal 
right-to-left exponentiation based on the binary representation of n, as the 
number of nonzero digits is on average only about one third of the total number 
of digits, but it is used here mainly for didactic reasons. What we take advan- 
tage here is the stunningly easy way to compute inverses in our group, after all, 
we only have to change the sign of the y-coordinate for this. 

Later on, we will also need the order of an elliptic curve, i.e., the number of its 
points. There are two essentially different ways of achieving this goal. The f irst 
one is to choose the coefficients a and b of the elliptic curve essentially at ran- 
dom and compute its order thereafter using Schoofs algorithm or its more ad- 
vanced version, the SEA-algorithm due to Schoof, Elkies and Atkins. We use 
here the second approach though, which makes use of a special class of elliptic 
curves, the so-called CM-curves. For this we need a set of discriminants for 
which the corresponding imaginary quadratic fields have a "small" class number. 
Below is a complete list of all those values belonging to the class numbers 1 or 2: 


D 
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A :z [-3, -4, -7, -8, -11, -19, -43, -67, -163, -15, -20, -24, -35, -40, 

-51, -52, -88, -91, -115, -123, -148, -187, -232, -235, -267, -403, -427] 

Next, we must choose for our given prime p the first number D from this set 
such that the equation 

4p = x 2 + | D |y 2 

has got an solution (x,y) with integers x and y. (In the unlikely case that no such 
D is available, we must switch to a different prime p.) In order to solve the Dio- 
phantine equation above for given values of p and D, we introduce another rou- 
tine due to Cornaccia-Smith: 

CS(p, D, a_, b_, c_, r_) := 

Prog 

If JACOBI (D, p) < 1 
RETURN false 
b_ := SQUARE_ROOT(D, p) 

If 0DD?(b_ - D) 
b_ :: p - b_ 
a_ 2*p 

c_ :z FLOOR(2.7p) 

Loop 

If b_ < c_ exit 
r_ :z MOD(a_, b_) 
a_ b_ 
b_ :z r_ 

a_ := 4-p - b_ A 2 
If MOD(a_, ABS(D)) > 0 
RETURN false 
c_ := J(a_/ABS(D)) 

If INTEGER?(c_) 

RETURN false 
[±b_, ±c_] 

[CS(23 , -3), CS(23 , -4), CS(23, -7)] 

[false, false, [+8, ±2]] 


I won't go into details as to how and why this routine works, but if you have a 
closer look at it, you will see, that it starts with an even square root of D mod p, 
(which might fail to exist though and leads to an abortion then) and contains 
elements of the Euclidean Algorithm thereafter which is carried out to certain 
point starting with 2p and the square root of D mod p above. The outcome might 
be used to construct the solutions (x,y) of the equation above, but there are 
still more chances of a failure, as you can see by inspecting the routine more 
closely. Even though, in view of the rather large list A the chances are good that 
we come up in the end with a handful of suitable values of D. 
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Now we must distinguish between three basic cases, namely 

I. D = -3 

II. D = -4 

III. D < -4 

We will assume the rather small value p = 37, so we might be able to check all 
computations manually. Furthermore, for this special value of p all 3 cases can be 
demonstrated. 

p 37 

SELECT (VECTOR?(CS(p , D)), D, A) 

[-3, -4, -7, -11, -67, -123] 


Case I: D=-3 


In the first place, there are six optional choices for the elliptic curves, which 
are all of the form 


y ! =x ! - g k . k=0,l,2,3,4,5 

for some quadratic nonresidue mod p. Moreover, in the case g = 1 mod 3 , which 
is true in this example, g must not be a cube mod p. For example, g = 2 is an ap- 
propriate value for p = 13 here due to 


(P - l)/3 

g := ITERATE(IFOAC0BI(g_, p) = -1 a MOD(g_ , p) U, g_, g_ + 1) , g_, 2) 

g := 2 

For each of these curves E it is true that #E is one of the six numbers 
p+l±u, p + l±(u± 3v)/2, where the signs of u and v may be chosen independ- 
ently here. In Order to find for each curve the correct value of n one must test 
the condition n-U = 0 for some randomly chosen point U t O on the elliptic curve. 
The following routine will do exactly this for us. 
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NOPl(p, k := 0, g_ := 1, n_, t_ p x_ := -l p u_ p v_, w_) := 

Prog 

w_ :z CS(p, -3) 

If - V ECTOR ?(w_) 

RETURN ? 

Loop 
g_ :+ 1 

If DACOBI(g_ p p) = -1 exit 
a :z 0 

b := HODSC- g_ A k p p) 
u_ := ABS( FIRST (w_)) 
v_ := ABS(FIRST(REST(w_))) 

t_ := {u_ p -u_ p (u_ + 3 ■ v_)/2 p (u_ - 3-v_)/2} 
t_ :z t_ u {- (u_ + 3-v_)/2 p - (u_ - 3-v_)/2} 

Loop 
Loop 
x_ : + 1 

If IACOBI(x_ a 3 + b p p) = 1 exit 
If x_ > p 

RETURN MAP_LIST(p + 1 + s_ p s_ p t_) 
w_ := [x_ p SQUARE_ROOT(x_ a 3 + b p p)] 

t_ := SELECT(mult(w_ p p + 1 + s_ p 0, p) = [p, p], s_ p t_) 
If DIM(t_) = 1 

RETURN p + 1 + FIRST(t_) 


We check this routine by computing the number of points on E: y 2 =x 3 -2 both by 
using it and a simple brute force computation: 


p - i 3 

NOP(a, b, p) :: p + 1 + I iACOBIfx + a*x + b, p) 

x=0 

r 2 3 k l 

VECTOR CLy = x + HODSC- 2 , p), n :z NOPl(p, k), n = NOP(a, b, p)J, k, 0, 5)' 

23 23 23 23 23 23 

y = x - ly = x - 2 y = x - 4 y = x - 8 y = x -16 y = x + 5 

n :z 48 n :z 49 n :z 39 n :z 28 n :z 27 n :z 37 

true true true true true true 


Case II: D=-4. 


Here the corresponding CM-curves are of the form 


y 2 =x 3 - g k x, k=0,l,2 ( 3 


where g is again some quadratic nonresidue mod p, for example g = 2 for p = 13 
as above. 
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This time, for each of these curves E the number #E of its points is given by 
one of the four numbers p + 1 ± u, p + 1 ± 2v, where (u,v) is again any solution of 
the Diophantine equation #4 above. The correct value among these can be found 
out in an analogous way as before using the following routine: 


NOP2(p, k := 0, g_ := 1, n_, t_, x_ := -1, u_, v_, w_) := 

Prog 

w_ := CS(p, -4) 

If - VECTOR? (w_) 

RETURN ? 

Loop 
g_ :+ 1 

If JACOBI(g_, p) = -1 exit 
a := HODS ( - g_ A k, p) 
b := 0 

u_ := ABS(FIRST(w_)) 
v_ := ABS(FIRST(REST(w_))) 
t_ := { u_ , — u_ , 2-v_, - 2-v_} 
x_ := 1 
Loop 
Loop 
x_ :+ 1 

If IACOBI(x_ a 3 + a-x_, p) = 1 exit 
If x_ > p 

RETURN MAP_LIST(p + 1 + s_, s_, t_) 
w_ := [x_, SQUARE_ROOT(x_ a 3 + a-x_, p)] 
t_ := SELECT(niult(w_, p + l + s_, a, p) = [p, p], s_, t_) 

If DIM(t_) = 1 

RETURN p + 1 + FIRST(t_) 

T 2 3 k 1 

VECTOR(Ly = x + MODS(- 2 , p)-x, n := NOP2(p, k), n = NOP(a, b, p)J, k, 0, 3)' 

"23 23 23 23 

y =x -x y =x - 2-x y = x - 4-x y = x - 8*x 

n := 40 n := 26 n := 36 n := 50 

true true true true 


Case III: D<-4. 


Here the corresponding CM-curves are of the form 
y 2 =x 3 -3rs 3 x+2rs 5 or y 2 =x 3 -3rs 3 g 2 x+2rs 5 g 3 


where g is again some quadratic nonresidue mod p, e.g. g = 2 for p = 13 as above. 
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But what about those mysterious numbers r and s occurring in these formulas, 
where do they come from? Well, telling you the whole story (just in case, you 
don't know it yet!) would take a lot of space here and lead us too far away. If you 
are interested in it (as well as in many other topics here I didn't have the time 
to dwell on) , I refer you to the wonderful book "Prime Numbers - A Compu- 
tational Perspective" by R. Crandall and C. Pomerance, which is sort of a "bible" 
when it comes to computational number theory. 

Here comes the solution: 



rslookup(D, p) :z 
Prog 

r := r 1 i st 4 POSITION ( D , rlist COL 1)|2 
r := r J.1 + r \2 .SQUARE_ROOT( r , p) 
s := s 1 i st 4 POSITION ( D , slist COL 1)|2 
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To simplify matters, as for those values of r and s, I'll simply look them up in the 
tables on the previous page, rather than compute them from scratch. Needless 
to say that these tables can also be found in the book quoted above though in a 
slightly different form, In particular, a triple [a,b,c] in the table named rlist 
corresponds to the term a + b e mod p, where again the square root must be 
computed mod p. 


For each of these curves E it is true that n = #E is one of the two numbers 
p + 1 ± u, where (u,v) is again a solution of our standard Diophantine equation #4 
above. The eventually correct value of n can be found in an analogous way as be- 
fore using the following routine: 

NOP3(p, D, k :: 0, g_ :z 1, n_, t_, x_ :: -1, u_, v_, w_) :z 
Prog 

w_ CS(p, D) 

If VECTOR? (w_) 

RETURN ? 
rslookupCD, p) 

Loop 
g_ :+ 1 

If IACOBI(g_, p) = -1 exit 
a := HODSC- 3 ■ r -s A 3 - g_ A (2 ■ k) , p) 
b := HODS(2 ■ r ■ s A 5 -g_ (3 ■ k) , p) 
u_ := ABS( FIRST (w_)) 
v_ := ABS( FIRST (REST (w_))) 
t_ := {u_ p -u_} 

Loop 
Loop 
x_ :+ 1 

If DACOBI(x_ a 3 + a*x_ + b, p) = 1 exit 
If x_ > p 

RETURN MAP_LIST(p + 1 + s_ p s_ p t_) 
w_ := [x_ p SQUAR E_ROOT ( x_ A 3 + a-x_ + b p p)] 
t_ := SELECT(mult(w_ p p + 1 + s_, a p p) = [p, p], s_ p t_) 

If DIM(t_) = 1 

RETURN p + 1 + FIRST (t_) 


REVERSECVECTORCLn := NOP3(p, -7, k) p n = NOP(a p b p p), y = x + a-x + bj p k, 0, 1)’) 

2 3 2 3 

y = x + 13 - x - 3 y = x + 15 * x + 13 

true true 

n := 32 n := 44 


Ok, one last routine and we are finally finished with our setup. It deals with the 
fact that not every elliptic curve is suitable for our purposes, but what we need 
is an elliptic curve of an appropriate order that is "almost" a prime, i.e., that is 
divisible only by small factors. The following routine will show if this condition is 
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sufficiently fulfilled by removing small factors from the order n of the curve up 
to a certain bound s. If it is a "good" curve, the remaining part should be a 
prime even if s is rather small. 

trialdiv(n, s :z 10 A 5, d_ :z 4, s_ :z [], t_ := 2, d_ :z 4) :z 
Loop 
Loop 

If MOD(n, t_) > 0 exit 
s_ := ADIOIN(t_, s_) 
n :/ t_ 

t_ :+ d_ A IF(t_ > 5) 
d_ :: 6 - d_ 

If t_ A 2 > ii 
If n > 1 

RETURN REVERSE(AD]OIN(ll, s_)) 

RETURN REVERSE(s_) 

If t_ > s 
Prog 

d_ := APPEND(IFCPRIHE?(n) p "P", "C"), STRING(DIM(n))) 

RETURN R EVERS E(AD JOIN (d_, s_)) 


Now, what do you think: Will our routines above pass the acid test when using 
primes of realistic order in a cryptographic environment, say with 160 bits and 
more? Let's try it out. 


159 159 

p := NEXT_PRIME(2 + RAN DOM (2 )) 

p := 979041373284628512978859729772372544317253699343 

SELECT (VECTOR? (CS(p , D)) p D p A) = [-3, -19, -163, -24, -123] 


Ok, we have in total 6+2+2+2+2=10 elliptic curves a tour disposal, let's check now 
if there is also one with a „big" prime factor. To be honest, there is also a lot of 
„garbage", but two or three are really good (see the accompanying DERIVE-File 
for a complete list!), the best one being the last below 


VECTOR(REVERSE(Ln := NOP3(p, -19, k) , trialdiv(n), y = x + a-x + bj), k, 0, 1) 
2 3 

y = x - 207368750592- x + 36382017816364032 [19, 467, C45] n := 

2 3 

. y = x - 1866318755328- x + 982314481041828864 [3, 3, 5, 11, 457, C43] n := 


979041373284628512978859508298659993048973669403 


979041373284628512978859951246085095585533729285 
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VECTORCREVERSECLn := N0P3(p, -123, k) , trialdiv(n) , y 


3 

X + 


,]>, 


a*x + bj) , k, 0, 1) 


2 3 

y = x + 108120918092691660115135537631025775915924068493 • x + 4319758727434414396267- 
2 3 

y = x - 5953110450403571942639891093140561073937082906* x - 85147915342623285823889- 
56571977728682405020466241 [11, 64763, P43] n := 

313869796106871491803609 [3, P48] n := 

979041373284628512978857852967500664853268633179 
979041373284628512978861606577244423781238765509 . 

As a matter of fact. I'll use exactly this one in part 2 of this treatise on the 
EIGamal dealing with all the goodies it has to offer. Hope to see you again here! 

Invitation to visit Europe’s Cultural Capital 2013 - Marseille 

Hello, 

I hope that you are keeping well . 

Maybe you will remember Bernard Eggger who was in Time Montreal in 2004 ? He is now president 
of the French association of math’s teachers for Provence Region (Aix-Marseille) 

And as Marseille is Europe’s capital for culture in 2013 , he decided to held the Association Annual 
Meeting there . 

I am also involved in the preparation of this meeting, not so much scientifically, but on the organisa- 
tion level. 

It is not exactly an international meeting because all the conferences and workshops will be in French 
,yet for this “capital” year it is planned to have a kind of Science Fair (called Souk) which will be 
more specifically open to teachers from abroad, if they which to present Posters, or whatever , even 
not in French. It is open to all of course, although we suppose that it will be more particularly of inter- 
est for teachers from Mediterranean countries. 

As it is the first (and probably last) time that I am busy myself with this kind of venture , 

I thought that I would let you know . 

So you may just visit the site : 

« Les mathematiques au carrefour des cultures de la Mediterranee » 

http://www.jnmarseille20 1 3 .fr 

With best regards 
Mit freundlichen GriiBen 
Marie-Faure Faurent 
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I prepared Piotr Trebisz’s 4 th part of his “Snail Shell” series. Guido Herweyers’ Statistics 4 
got its self dynamic and needed more space than expected. So I must leave the Conical 
Shells for the next DNL. See some “mouthwatering” plots. 
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Impressions from former Indochina 
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Impressions from Laos, Vietnam and Cambodia 



Sunset over the Mekong, Luang Prabang, Laos 


Vientiane, Laos 





One Pillar Pagoda, Hanoi 


Students in the Temple of Literature, Hanoi 



Halong Bay 
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Bay on, Angkor Wat, Cambodia 


Calculus in Cambodia 





